#define _CRT_SECURE_NO_DEPRECATE
#include <cstdio>
#include <fcntl.h> // for _O_BINARY
#include <vector>
#include <string>
#include <cstring>
#include <algorithm>
#include <cmath>
#include "grate_constants.h"
#include "ext_rational.h"
#include "words_queue.h"

using namespace std;

// Parameters that are read from the ini-file, or set automatically, or defined here (values by default)
bool auto_queue_size = true;
bool auto_dfa_size = true;

int QUEUE_SIZE;
int DFA_SIZE; // less than 2^28 !!!!
int MAXITER = 20000; // Iterations the algorithm with oracle does

bool USE_ORACLE = false; // true - Algo1Q, false - Algo1
bool ADD_LOOPS = false; // if we add loops, it makes algorithm work slower, but DFA is always primitive
bool IMPROVE_TARJAN = false; // if true, assumes that the resulting dfa always have one scc, this saves memory
int GEN_EDGE_LIST = -1; // replaces DFA with edge list, sometimes this works faster (0 - don't use, 1 - use, -1 - auto)

bool OUTPUT_PROGRESS = true; // writes additional data to the log file

// Input data:
int ALPH_SIZE;
int ROOT_LENGTH;
ext_rational beta;
enum {WEAK, MODERATE, STRONG} power_type; // 0 - weak, 1 - moderate, 2 - strong
double delta;

FILE * logfile;

// Main containers - queue and DFA table
words_queue * WORDS_QUEUE;
// todo
int total;
int * dfa;

//------------------------------------------- INITIALIZATION -----------------------------------

// Determines if queue size is more important than dfa size
int queue_more_important(){
  if (ALPH_SIZE <= 2) return 0;
  return (beta < ext_rational(ALPH_SIZE-1, ALPH_SIZE-2, 1));
}

void create_queue(int size){
  int qsize;
  if (size != -1){
    qsize = size;
  }
  else{ // define size automatically
    double max_size = ROOT_LENGTH * pow((double)ALPH_SIZE, (double)ROOT_LENGTH);
    if (sizeof(char) * max_size < LITTLE_MEMORY){
      qsize = (int)max_size + 1;
    }
    else{
      qsize = QUEUE_MEMORY_STEP;
      while (1){
        try{
          char * buff = new char [qsize + QUEUE_MEMORY_STEP];
          delete[] buff;
          qsize += QUEUE_MEMORY_STEP;
        }
        catch(...){
          break;
        }
      }
      qsize >>= 1;
    }
  }
  //words_queue = new char [QUEUE_SIZE];
  bool use_long_mask = (ALPH_SIZE > 4);
  WORDS_QUEUE = new words_queue(qsize, use_long_mask);
}

void set_dfa_size(){
  double max_size = ROOT_LENGTH * pow((double)ALPH_SIZE, (double)ROOT_LENGTH);
  if (sizeof(int) * max_size * ALPH_SIZE < LITTLE_MEMORY){
    DFA_SIZE = (int)max_size + 1;
  }
  else{
    DFA_SIZE = DFA_MEMORY_STEP;
    while (1){
      try{
        dfa = new int [(DFA_SIZE + DFA_MEMORY_STEP) * ALPH_SIZE];
        delete[] dfa;
        DFA_SIZE += DFA_MEMORY_STEP;
      }
      catch(...){
        break;
      }
    }
    DFA_SIZE >>= 1;
  }
  dfa = new int [DFA_SIZE * ALPH_SIZE];
}

void init(){
  if (auto_queue_size == false){
    create_queue(QUEUE_SIZE);
  }
  if (auto_dfa_size == false){
    dfa = new int [DFA_SIZE * ALPH_SIZE];
  }
  
  if (queue_more_important()){
    if (auto_queue_size) create_queue(-1);
    if (auto_dfa_size) set_dfa_size();
  }
  else{
    if (auto_dfa_size) set_dfa_size();
    if (auto_queue_size) create_queue(-1);
  }

  if (OUTPUT_PROGRESS){
    if (auto_queue_size) fprintf(logfile, "Amount of memory allocated for queue = %d bytes\n", sizeof(char)*WORDS_QUEUE->size);
    if (auto_dfa_size) fprintf(logfile, "Amount of memory allocated for DFA = %d bytes\n", sizeof(int)*DFA_SIZE*ALPH_SIZE);
  }

  total = 1;
  memset(dfa, 0, sizeof(int) * DFA_SIZE * ALPH_SIZE);
}


// ------------------------------ ANTIDICTIONARY CONSTRUCTION --------------------------

i64 max_words_queue_size;
int words_num; // Number of words in the DFA

// Returns the miminal length of forbidden word where plen is the length of its root
inline int get_length(int plen){
  return (2 * plen * beta.num + beta.p + 2 * beta.denum - 1) / (2 * beta.denum);
}

// Adds the word to automaton
void add_word(const char * word){
  if (OUTPUT_ANTIDICTIONARY){
    printf("%s\n", word);
  }
  int cv = 1, i;
  words_num++;
  for (i=0; word[i+1]; i++){
    int c = word[i]-'a';
    if (!dfa[cv*ALPH_SIZE + c]){
      total++;
      if (total >= DFA_SIZE){
        fprintf(logfile, "ERROR! Insufficient DFA size!\n");
        exit(1);
      }
      dfa[cv*ALPH_SIZE + c] = total;
    }
    cv = dfa[cv*ALPH_SIZE + c];
  }
  dfa[cv*ALPH_SIZE + word[i]-'a'] = -1;
}

// Checks if some suffix of the word in buff is forbidden
bool bad_suffix_trie(const char * buff, int len){
  int cv = 1, i, cur = 0;
  char bij[MAX_ALPH_SIZE];
  memset(bij, 0, sizeof(char)*ALPH_SIZE);
  for (i=len-1; i>=0; i--){
    if (!bij[buff[i]-'a']){
      bij[buff[i]-'a'] = 'a' + cur;
      cur++;
    }
    int c = bij[buff[i]-'a'] - 'a';
    int zv = dfa[cv * ALPH_SIZE + c];
    if (!zv) return false;
    if (zv < 0) return true;
    cv = zv;
  }
  return false;
}

bool check_head_tail(const char * buff, int len, int head, int tail, int zquant_letters[MAX_WORD_LENGTH][MAX_ALPH_SIZE]){
  int hlen = tail - head;
  int tlen = len - tail;
  if (hlen * beta.intval() + tlen > len) return false;
  int tail_quant[MAX_ALPH_SIZE];
  int head_quant[MAX_ALPH_SIZE];
  int qmax[MAX_ALPH_SIZE];
  for (int j=0; j<ALPH_SIZE; j++){
    tail_quant[j] = zquant_letters[len][j] - zquant_letters[tail][j];
    head_quant[j] = zquant_letters[head + hlen][j] - zquant_letters[head][j];
    if (head_quant[j] < tail_quant[j]) return false;
    qmax[j] = zquant_letters[head + tlen][j] - zquant_letters[head][j];
    if (beta.intval()== 1 && power_type == STRONG){
      if (qmax[j] != tail_quant[j]) return false;
    }
  }
  for (int i=1; i<beta.intval(); i++){
    for (int j=0; j<ALPH_SIZE; j++){
      int cquant;
      // Each part is a perm of "head"
      cquant = zquant_letters[head - i * hlen + hlen][j] - zquant_letters[head - i * hlen][j];
      if (cquant != head_quant[j]) return false;
      cquant = zquant_letters[head - i * hlen + tlen][j] - zquant_letters[head - i * hlen][j];
      if (cquant > qmax[j]) qmax[j] = cquant;
      if (i == beta.intval() - 1 && power_type == STRONG){
        if (cquant != tail_quant[j]) return false;
      }
    }
  }
  if (power_type == MODERATE){
    for (int j=0; j<ALPH_SIZE; j++){
      if (qmax[j] < tail_quant[j]) return false;
    }
  }
  return true;
}

bool bad_suffix_check(const char * buff, int len, int zquant_letters[MAX_WORD_LENGTH][MAX_ALPH_SIZE]){
  int cnt[MAX_ALPH_SIZE];
  memset(cnt, 0, sizeof(int) * MAX_ALPH_SIZE);
  //calc_zzquant(buff, len);
  int tail = len, head = len-1;
  while (head >= 0){
    tail--;
    char c = buff[tail] - 'a';
    if (tail == len-1) cnt[c]--;
    else cnt[c] -= 2;
    while (cnt[c] < 0){
      head--;
      if (head < 0) break;
      cnt[buff[head] - 'a']++;
    }
    if (head >= 0){
      int zhead = head;
      while (zhead >= 0){
        int hlen = tail - zhead;
        int tlen = len - tail;
        if (get_length(hlen) > hlen * beta.intval() + tlen) break;
        if (power_type == WEAK && beta.intval() == 1) return true;
        // ZQUANT_LETTERS GOES INTO WORK HERE!!!!
        if (check_head_tail(buff, len, zhead, tail, zquant_letters)) return true;
        zhead--;
      }
    }
  }
  return false;
}

// Checks if there is a bad suffix in the word in buff
// Calls check in a trie and then kind of dumb check
bool bad_suffix(const char * buff, int len, int zquant_letters[MAX_WORD_LENGTH][MAX_ALPH_SIZE]){
  if (beta.num % beta.denum == 0 && beta.p == 0){
    return bad_suffix_trie(buff, len);
  }
  else if (beta.intval() == 1 && power_type == STRONG){
    return bad_suffix_trie(buff, len) || bad_suffix_check(buff, len, zquant_letters);
  }
  else{
    return bad_suffix_check(buff, len, zquant_letters);
  }
}

int root_length;
char root_letnum;
char root_letters[MAX_ALPH_SIZE];
int root_freqs[MAX_ALPH_SIZE];
char cur_letters[MAX_ALPH_SIZE];
int cur_freqs[MAX_ALPH_SIZE];

// quant_letters[i][j] - how many letters "j" we have in the prefix of length "i";
int quant_letters[MAX_WORD_LENGTH][MAX_ALPH_SIZE];

void build_words(char* cword, int clen, int nlen, int mlen, char letnum){
  if (clen >= nlen && clen <= mlen){
    bool should_add;
    if (power_type == WEAK){
      should_add = true;
    }
    else if (power_type == STRONG){
      int mod = clen % root_length;
      int i;
      // check that the last "mod" symbols are permutations of the first "mod"
      for (i=0; i<ALPH_SIZE; i++){
        if (quant_letters[clen][i] - quant_letters[clen-mod][i] != quant_letters[mod][i]) break;
      }
      should_add = (i == ALPH_SIZE);
    }
    else if (power_type == MODERATE){
      int mod = clen % root_length;
      int i;
      // check that in the last "mod" we have not greater symbols of each kind
      // than in the first "mod" symbols of each part
      for (i=0; i<ALPH_SIZE; i++){
        int hmuch = quant_letters[clen][i] - quant_letters[clen-mod][i];
        int maxq = 0;
        for (int j=0; j<beta.intval(); j++){
          int curq = quant_letters[root_length * j + mod][i] - quant_letters[root_length * j][i];
          if (curq > maxq) maxq = curq;
        }
        if (hmuch > maxq) break;
      }
      should_add = (i == ALPH_SIZE);
    }
    if (should_add){
      add_word(cword);
      return;
    }
  }
  if (clen == mlen || bad_suffix(cword, clen, quant_letters)){
    return;
  }
  char t_letters[MAX_ALPH_SIZE];
  int t_freqs[MAX_ALPH_SIZE];
  if (clen % root_length == 0){
    letnum = root_letnum;
    for (int i=0; i<ALPH_SIZE; i++){
      t_letters[i] = cur_letters[i];
      cur_letters[i] = root_letters[i];
      t_freqs[i] = cur_freqs[i];
      cur_freqs[i] = root_freqs[i];
    }
  }
  for (char c=0; c<letnum; c++){
    cword[clen] = cur_letters[c];
    cword[clen + 1] = 0;
    bool ok = false;
    for (char j=0; j<ALPH_SIZE; j++){
      quant_letters[clen + 1][j] = quant_letters[clen][j];
    }
    quant_letters[clen + 1][cword[clen] - 'a']++;
    if (clen < nlen - 1) ok = !bad_suffix(cword, clen + 1, quant_letters);
    else ok = !bad_suffix(cword + 1, clen, &(quant_letters[1]));
    if (ok){
      cur_freqs[c]--;
      bool swapped = false;
      if (cur_freqs[c] == 0){
        swapped = true;
        swap(cur_freqs[c], cur_freqs[letnum-1]);
        swap(cur_letters[c], cur_letters[letnum-1]);
        build_words(cword, clen + 1, nlen, mlen, letnum - 1);
      }
      else{
        build_words(cword, clen + 1, nlen, mlen, letnum);
      }
      if (swapped){
        swap(cur_freqs[c], cur_freqs[letnum-1]);
        swap(cur_letters[c], cur_letters[letnum-1]);
      }
      cur_freqs[c]++;
    }
    for (char j=0; j<ALPH_SIZE; j++){
      quant_letters[clen + 1][j] = 0;
    }
  }
  cword[clen] = 0;
  if (clen % root_length == 0){
    for (int i=0; i<ALPH_SIZE; i++){
      cur_letters[i] = t_letters[i];
      cur_freqs[i] = t_freqs[i];
    }
  }
}

// Hadles root stored in cword and adds forbidden words constructed from it if necessary
void handle_word(const char* cword, int clen, char maxlet){
  int i;
  char thword[MAX_WORD_LENGTH + 1];

  root_length = clen;
  root_letnum = maxlet-'a'+1;

  memset(root_freqs, 0, sizeof(int) * MAX_ALPH_SIZE);

  for (char ch='a'; ch<=maxlet; ch++){
    root_letters[ch-'a'] = ch;
  }

  // Copy the root to thword
  for (i=0; i<root_length; i++){
    thword[i] = cword[i];
    root_freqs[cword[i]-'a']++;
  }
  thword[clen] = 0;
  if (clen > MAX_WORD_LENGTH){
    fprintf(logfile, "Achtung! The length of the word is too small!");
    exit(1);
  }

  // The length we need
  int nlen = get_length(clen);
  // Maximal possible length of forbidden word for strict powers
  int maxlen;
  if (power_type == WEAK){
    maxlen = nlen;
  }
  else{
    maxlen = (nlen + root_length - 1) / root_length * root_length;
  }

  build_words(thword, clen, nlen, maxlen, root_letnum);  
}


// Main trie build procedure
void build_antidictionary(){
  long long marker_pos;
  
  WORDS_QUEUE->push_marker(marker_pos); // sets marker_pos to 0
  WORDS_QUEUE->push_symbol('a');

  int i;
  int clen = 0;
  char bad_word[MAX_WORD_LENGTH + 1];

  while (! WORDS_QUEUE->is_empty()){
    if (WORDS_QUEUE->is_index(marker_pos)){
      if (clen < ROOT_LENGTH){
        clen++;
        WORDS_QUEUE->push_marker(marker_pos);
        WORDS_QUEUE->pop_symbol();
        continue;
      }
      else break;
    }
    
    // there is a word in the queue starting at qbeg with (length = clen)
    char maxlet = 'a';
    for (int j=0; j<ALPH_SIZE; j++){
      quant_letters[0][j] = 0;
    }
    for (i=0; i<clen; i++){
      bad_word[i] = WORDS_QUEUE->pop_symbol() + 'a';
      if (bad_word[i] > maxlet){
        maxlet = bad_word[i];
      }
      for (int j=0; j<ALPH_SIZE; j++){
        quant_letters[i+1][j] = quant_letters[i][j];
      }
      quant_letters[i+1][bad_word[i]-'a']++;
    }
    bad_word[clen] = 0;

    handle_word(bad_word, clen, maxlet);

    if (clen < ROOT_LENGTH){
      char ch;
      for (ch='a'; ch<=maxlet; ch++){
        bad_word[clen] = ch;
        bad_word[clen+1] = 0;
        for (int j=0; j<ALPH_SIZE; j++){
          quant_letters[clen+1][j] = quant_letters[clen][j];
        }
        quant_letters[clen+1][ch-'a']++;
        if (!bad_suffix(bad_word, clen+1, quant_letters)){
          if (!WORDS_QUEUE->push_word(bad_word)){
            fprintf(logfile, "ERROR! Insufficient words queue size!\n");
            exit(1);
          }
        }
      }
      if (maxlet-'a'+1 < ALPH_SIZE){
        bad_word[clen] = maxlet+1;
        bad_word[clen+1] = 0;
        if (!WORDS_QUEUE->push_word(bad_word)){
            fprintf(logfile, "ERROR! Insufficient words queue size!\n");
            exit(1);
          }
      }
    }
    i64 cwqsz = WORDS_QUEUE->get_number_of_letters();
    if (cwqsz > max_words_queue_size){
      max_words_queue_size = cwqsz;
    }
  }

  delete WORDS_QUEUE;
}



// ---------------------------- AUTOMATON CONSTRUCTION ----------------------------

//char * last;
//int * lf;
const int LAST_MASK = (15 << 28);
const int LF_MASK = (1<<28) - 1;
int * lf_last;

char * sigma;

int back_pack_mode; // 0 - 4 bits per mask, 1 - 8 bits per mask, 2 - 16 bits per mask
short * back;

int * vert_queue;
int vqbeg, vqend;
int vqsize; // always a perfect power of 2

inline void push_vertex(int x){
  vert_queue[vqend] = x;
  vqend = (vqend + 1) & (vqsize - 1);
}

int SIGMA_LEN;

inline char get_sval(char * s, int c){
  return (c & 1) ? ((s[c>>1] & (LONG_MASK << 4)) >> 4) : (s[(c>>1)] & LONG_MASK);
}

inline void set_sval(char * s, int c, char d){
  if (c & 1){
    s[c >> 1] &= LONG_MASK;
    s[c >> 1] |= (d << 4);
  }
  else{
    s[c >> 1] &= (LONG_MASK << 4);
    s[c >> 1] |= d;
  }
}

void make_permut(char* p1, char* p2){
  int i;
  for (i=0; i<ALPH_SIZE; i++){
    if (get_sval(p1, i) != LONG_MASK){
      set_sval(p1, i, get_sval(p2, get_sval(p1, i)));
    }
  }
}

inline char get_sigma(int v, int c){
  return get_sval(sigma + v*SIGMA_LEN, c);
}

inline void set_sigma(int v, int c, char d){
  set_sval(sigma + v*SIGMA_LEN, c, d);
}

inline void set_last(int v, char c){
  lf_last[v] &= ~LAST_MASK;
  lf_last[v] |= ((int)(c) << 28);
}

inline void set_lf(int v, int clf){
  lf_last[v] &= ~LF_MASK;
  lf_last[v] |= clf;
}

inline char get_last(int v){
  char z = (((lf_last[v] & LAST_MASK) >> 28) & LONG_MASK);
  return (z == 15) ? -1 : z;
}

inline int get_lf(int v){
  return lf_last[v] & LF_MASK;
}

inline void set_back(int v, int c){
  if (back_pack_mode == 0){
    int mod = v&3;
    back[v>>2] |= (1 << (4*mod + c));
  }
  else if (back_pack_mode == 1){
    int mod = v&1;
    back[v>>1] |= (1 << (8*mod + c));
  }
  else{
    back[v] |= (1 << c);
  }
}

inline bool is_back(int v, int c){
  if (back_pack_mode == 0){
    int mod = v&3;
    return (back[v>>2] & (1 << (4*mod + c))) > 0;
  }
  else if (back_pack_mode == 1){
    int mod = v&1;
    return (back[v>>1] & (1 << (8*mod + c))) > 0;
  }
  else{
    return (back[v]  & (1 << c)) > 0;
  }
}

void cnt_vertex(int v){
  int c;
  char * zsigma = new char [SIGMA_LEN];
  char clast = get_last(v);
  for (c=0; c<=clast+1; c++){
    if (c >= ALPH_SIZE) break;
    if (!dfa[v*ALPH_SIZE + c]){
      char d = get_sigma(v, c);
      if (d == LONG_MASK){
        d = get_last(get_lf(v)) + 1;
      }
      dfa[v*ALPH_SIZE + c] = dfa[get_lf(v)*ALPH_SIZE + d];
      set_back(v, c);
    }
    else{
      int u = dfa[v*ALPH_SIZE + c];
      if (u < 0){
        dfa[v*ALPH_SIZE + c] = 0;
        continue;
      }
      push_vertex(u);
      set_last(u, (c > clast) ? c : clast);
      memcpy(sigma + u*SIGMA_LEN, sigma + v*SIGMA_LEN, sizeof(char)*SIGMA_LEN);
      char d = get_sigma(v, c);
      if (d == LONG_MASK){
        d = get_last(get_lf(v)) + 1;
        set_sigma(u, c, d);
      }
      set_lf(u, dfa[get_lf(v)*ALPH_SIZE + d]);
      int zv = get_lf(v);
      while (is_back(zv, get_sigma(u, c))){
        memcpy(zsigma, sigma + zv*SIGMA_LEN, sizeof(char) * SIGMA_LEN);
        if (get_sval(zsigma, get_sigma(u, c)) == LONG_MASK){
          set_sval(zsigma, get_sigma(u, c), get_last(get_lf(zv)) + 1);
        }
        make_permut(sigma + u*SIGMA_LEN, zsigma);
        zv = get_lf(zv);
      }
    }
  }
  for (c=clast+2; c<ALPH_SIZE; c++){
    dfa[v*ALPH_SIZE + c] = dfa[v*ALPH_SIZE + clast+1];
    if (is_back(v, clast+1)){
      set_back(v, c);
    }
  }
  delete[] zsigma;
}

// !!!! --------------   DEBUGGING PURPOSES!!! --------------------------
/*void read_word(const char * word){
  int cv = 1, i;
  for (i=0; word[i]; i++){
    int c = get_sigma(cv, word[i] - 'a');
    if (!dfa[cv*ALPH_SIZE + c]){
      c=c;
    }
    else cv = dfa[cv*ALPH_SIZE + c];
  }
}*/

void build_dfa(){
  
  SIGMA_LEN = (ALPH_SIZE + 1) / 2;
  sigma = new char [(total + 1) * SIGMA_LEN];

  int BACK_SIZE;

  if (ALPH_SIZE <= 4){
    back_pack_mode = 0;
    BACK_SIZE = (total + 4) / 4;
  }
  else if (ALPH_SIZE <= 8){
    back_pack_mode = 1;
    BACK_SIZE = (total + 2) / 2;
  }
  else{
    back_pack_mode = 2;
    BACK_SIZE = (total + 1);
  }
  back = new short [BACK_SIZE];
  memset(back, 0, sizeof(short) * BACK_SIZE);

  lf_last = new int [total + 1];
  int i;
  for (i=1; i<=total; i++){
    lf_last[i] = 0;
  }
  memset(sigma, 0xFF, sizeof(char) * (total + 1) * SIGMA_LEN);
  for (i=1; i<ALPH_SIZE; i++){
    dfa[1*ALPH_SIZE + i] = 2;
  }
  
  vqsize = 1;
  while (vqsize < words_num) vqsize <<= 1;
  vert_queue = new int [vqsize];
  vqbeg = 0; vqend = 0;
  
  push_vertex(2);
  set_last(1, -1);
  set_lf(1, 1);
  set_lf(2, 1);

  while (vqbeg != vqend){
    cnt_vertex(vert_queue[vqbeg]);
    vqbeg = (vqbeg + 1) & (vqsize - 1);
  }

  // DEBUGGING PURPOSES !!!!
  //read_word("aaabababbbabaaababbbaaabbbaabbbabaaab");

  delete[] sigma;
  delete[] back;
  delete[] vert_queue;
  delete[] lf_last;
}

//--------------------------------------------------- ANSWER CALCULATION-----------------------------------------

double * cnt_new;
double * cnt_old;

double grate;

void cnt_algo2(){
  int i, j, k;
  char h;
  cnt_old = new double [total+1];
  memset(cnt_old, 0, sizeof(double) * (total + 1));
  if (OUTPUT_PROGRESS){
    fprintf(logfile, "old counters allocated\n");
    fflush(logfile);
  }
  cnt_new = new double [total+1];
  memset(cnt_new, 0, sizeof(double) * (total + 1));
  if (OUTPUT_PROGRESS){
    fprintf(logfile, "new counters allocated\n");
    fflush(logfile);
  }
  cnt_old[1] = 1;
  double pcomp = 1, ccomp, last_loc_max = 0, last_loc_min = 0;
  double * rates = new double [MAXITER];
  for (i=0; i<MAXITER; i++){
    for (j=1; j<=total; j++){
      for (h=0; h<ALPH_SIZE; h++){
        k = dfa[j*ALPH_SIZE + h];
        if (!k) continue;
        cnt_new[k] += cnt_old[j];
      }
    }
    ccomp = 0;
    for (j=1; j<=total; j++){
      if (ADD_LOOPS) cnt_new[j] += cnt_old[j];
      ccomp += cnt_new[j];
      cnt_old[j] = cnt_new[j];
      cnt_new[j] = 0;
    }
    rates[i] = ccomp / pcomp;
    // ORACLE HERE
    if (i >= MINITER){
      if ((rates[i-1] >= rates[i-2]) && (rates[i-1] >= rates[i])) last_loc_max = rates[i-1];
      if ((rates[i-1] <= rates[i-2]) && (rates[i-1] <= rates[i])) last_loc_min = rates[i-1];
      if (fabs(rates[i] - rates[i-1]) < delta){
        // Oscillating
        if ((last_loc_max > 0) && (last_loc_min > 0)){
          if (last_loc_max - last_loc_min < delta){
            if (OUTPUT_PROGRESS){
              fprintf(logfile, "Oscillating sequence! The number of iterations = %d\n", i);
            }
            if (GUI_VERSION) printf("%d\n", i);
            grate = rates[i] - (ADD_LOOPS == true);
            return;
          }
        }
        // Monotonic
        else{
          if (OUTPUT_PROGRESS){
            fprintf(logfile, "Monotonic sequence! The number of iterations = %d\n", i);
          }
          if (GUI_VERSION) printf("%d\n", i);
          double q = pow(fabs(rates[i]-rates[i-1]), 1.0 / i);
          double sum = (rates[i] - rates[i-1]) / (1 - q);
          grate = rates[i-1] + sum - (ADD_LOOPS == true);
          return;
        }
      }
    }
    pcomp = ccomp;
  }
  if (OUTPUT_PROGRESS){
    fprintf(logfile, "\nNot enough iterations! Current difference is %.20lf\n\n", rates[i-1] - rates[i-2]);
  }
  if (GUI_VERSION) printf("-1\n");
  grate = rates[i-1] - (ADD_LOOPS == true);
}

// Edge list
int * edge1;
int * edge2;
bool have_edge_list;

// Some values to output
int scc_num, nontriv_scc_num, scc_size, scc_edgnum;

int * zcc;

void cnt_algo1(int id){
  double minrate, maxrate, crate, alpha;
  int i, iter = 0;
  while (iter < MAXITER){
    iter++;
    if (have_edge_list){
      for (i=0; i<scc_edgnum; i++){
        cnt_new[edge2[i]] += cnt_old[edge1[i]];
      }
    }
    else{
      /*for (i=1; i<=scc_size; i++){
        for (j=0; j<ALPH_SIZE; j++){
          if (dfa[i*ALPH_SIZE + j] > 0){
            cnt_new[dfa[i*ALPH_SIZE + j]] += cnt_old[i];
          }
        }
      }*/
    }
    minrate = 1e100;
    maxrate = 0;
    for (i=1; i<=scc_size; i++){
      if (ADD_LOOPS){
        cnt_new[i] += cnt_old[i];
      }
      crate = cnt_new[i] / cnt_old[i];
      if (crate < minrate) minrate = crate;
      if (crate > maxrate) maxrate = crate;
      cnt_old[i] = cnt_new[i];
      cnt_new[i] = 0;
    }
    if (maxrate - minrate < 2 * delta){
      alpha = (minrate + maxrate) / 2 - (ADD_LOOPS == true);
      if (minrate > MAX_GRATE) alpha = 0;
      if (alpha > grate) grate = alpha;
      if (OUTPUT_PROGRESS) fprintf(logfile, "\n%d iterations used\n\n", iter);
      fprintf(logfile, "Growth rate at scc %d: %.20lf\n", id, alpha);
      if (GUI_VERSION) printf("%d\n", iter);
      return;
    }
  }
  if (OUTPUT_PROGRESS) fprintf(logfile, "\nNot enough iterations or precision error! %.20lf precision only!\n\n", (maxrate-minrate)/2);
  if (GUI_VERSION) printf("-1\n");
  alpha = (maxrate + minrate) / 2 - (ADD_LOOPS == true);
  fprintf(logfile, "Growth rate at scc %d: %.20lf\n", id, alpha);
  if (alpha > grate) grate = alpha;
  
}

void solve_scc(int id){
  int i, j;

  scc_edgnum = 0;
  scc_size = 0;
  
  for (i=1; i<=total; i++){
    if (zcc[i] == id){
      scc_size++;
      zcc[i] = -scc_size;
    }
  }

  for (i=1; i<=total; i++){
    if (zcc[i] < 0){
      for (j=0; j<ALPH_SIZE; j++){
        int v = dfa[i*ALPH_SIZE + j];
        if (v > 0 && zcc[v] < 0){
          scc_edgnum++;
        }
      }
    }
  }

  have_edge_list = true;
  
  edge1 = new int [scc_edgnum];
  edge2 = new int [scc_edgnum];
  int cnum = 0, znum = 0;
  for (i=1; i<=total; i++){
    if (zcc[i] < 0){
      for (j=0; j<ALPH_SIZE; j++){
        int v = dfa[i*ALPH_SIZE + j];
        if (v > 0 && zcc[v] < 0){
          edge1[cnum] = -zcc[i];
          edge2[cnum++] = -zcc[v];
        }
      }
    }
  }

  for (i = 1; i <= total; i++){
    if (zcc[i] < 0) zcc[i] = 0;
  }

  cnt_old = new double [scc_size + 1];
  if (OUTPUT_PROGRESS){
    fprintf(logfile, "old counters allocated\n");
    fflush(logfile);
  }
  cnt_new = new double [scc_size + 1];
  if (OUTPUT_PROGRESS){
    fprintf(logfile, "new counters allocated\n");
    fflush(logfile);
  }
  
  for (i=1; i<=scc_size; i++){
    cnt_new[i] = 0;
    cnt_old[i] = 1;
  }
  
  if (scc_edgnum == scc_size){
    if (grate == 0) grate = 1;
    if (have_edge_list){
      delete[] edge1;
      delete[] edge2;
    }
    return;
  }
  
  nontriv_scc_num++;
  cnt_algo1(id);
  if (have_edge_list){
    delete[] edge1;
    delete[] edge2;
  }
}

// TARJAN ALGORITHM

int * mtime;
int * father;
int * stack;
int ssize;
char * nxt_used; // binary encodeing here: 15-mask of number of edges done [+16] = used
char * scc;

void dfs(){
  memset(nxt_used, 0, sizeof(char) * (total + 1));
  memset(mtime, 0, sizeof(int) * (total + 1));
  ssize = 0;
  father[1] = 0;
  mtime[1] = 1;
  int ctime = 1, cver = 1, cedge, nver, pver;
  while (cver > 0){
    cedge = nxt_used[cver] & LONG_MASK;
    if (cedge == 0){ 
      nxt_used[cver] ^= 16; // is in stack
      stack[++ssize] = cver;
    }
    else{ // check the son
      pver = dfa[cver*ALPH_SIZE + cedge-1];
      if ((pver > 0) && (father[pver] == cver) && (mtime[pver] < mtime[cver])){
        mtime[cver] = mtime[pver];
        nxt_used[cver] |= 32; // improvement
      }
    }
    if (cedge >= ALPH_SIZE){
      // handle vertex away

      if ((nxt_used[cver] & 32) == 0){ // no improvement
        if (stack[ssize] == cver){
          nxt_used[cver] ^= 16; // ouf of stack
          ssize--;
        }
        else{
          scc_num++;
          do{
            pver = stack[ssize--];
            nxt_used[pver] ^= 64; // in scc
            scc_size++;
            nxt_used[pver] ^= 16; // ouf of stack
            scc[pver] = scc_num;
          }while (pver != cver);
        }
      }

      cver = father[cver];
      continue;
    }
    nxt_used[cver]++;
    nver = dfa[cver*ALPH_SIZE + cedge];
    if (!nver) continue;
    if (!mtime[nver]){
      ctime++;
      mtime[nver] = ctime;
      father[nver] = cver;
      cver = nver;
      continue;
    }
    else if ((nxt_used[nver] & 16) && (mtime[nver] < mtime[cver])){
      mtime[cver] = mtime[nver];
      nxt_used[cver] |= 32; // improvement
    }
  }
}

void divide_scc_solve(){
  if (!USE_ORACLE){
    mtime = new int [total + 1];
    father = new int [total + 1];
    stack = new int [total + 1];
    nxt_used = new char [total + 1];
    
    scc = new char [total + 1];
    memset(scc, 0, sizeof(char) * (total + 1));
    
    dfs();
    
    delete[] mtime;
    delete[] stack;
    delete[] father;
    delete[] nxt_used;
    
    if (OUTPUT_PROGRESS){
      fprintf(logfile, "The number of non-singleton strongly connected components: %d\n", scc_num);
      fflush(logfile);
    }

    zcc = new int[total + 1];
    for (int i = 1; i <= total; i++){
      zcc[i] = scc[i];
    }
    delete[] scc;

    for (int cur_scc=1; cur_scc<=scc_num; cur_scc++){
      solve_scc(cur_scc);
      if (OUTPUT_PROGRESS){
        fprintf(logfile, "The number of vertices in scc %d: %d\n", cur_scc, scc_size);
        fprintf(logfile, "The number of edges in scc %d: %d\n", cur_scc, scc_edgnum);
        fprintf(logfile, "\n");
        fflush(logfile);  
      }
    }

    delete[] zcc;

    if (OUTPUT_PROGRESS){
      fprintf(logfile, "The number of non-trivial scc: %d\n", nontriv_scc_num);
      fflush(logfile);
    }

  }
  else{
    cnt_algo2();
  }
}

// ------------------------------------ CLEAR DATA -------------------------

void clear_data(){
  grate = 0;
  scc_num = 0;
  nontriv_scc_num = 0;
  scc_size = 0;
  scc_edgnum = 0;
  words_num = 0;
  max_words_queue_size = 0;
  
  /*if (USE_ORACLE){
    delete[] dfa;
  }
  if (!USE_ORACLE && !have_edge_list){
    delete[] dfa;
  }*/
  delete[] dfa;

  delete[] cnt_new;
  delete[] cnt_old;
}

void solve(){

  char logname[100];
  char sgn = (beta.p == 0)?'-':'+';
  sprintf(logname, "a%dr%db%d%d%c%s.log", ALPH_SIZE, ROOT_LENGTH, beta.num, beta.denum, sgn, 
    (power_type == 0) ? "weak" : (power_type == 1) ? "moder" : "strong");
  logfile = fopen(logname, "w");

  if (OUTPUT_PROGRESS){
    fprintf(logfile, "Calculating answer for test:\n");
    fprintf(logfile, "Alphabet size = %d, Root length <= %d\nExponent = %d/%d%c", ALPH_SIZE, ROOT_LENGTH, beta.num, beta.denum, sgn);
    char type[20];
    if (power_type == WEAK) strcpy(type, "WEAK");
    else if (power_type == MODERATE) strcpy(type, "MODERATE");
    else strcpy(type, "STRONG");
    fprintf(logfile, "%s", type);
    fprintf(logfile, "\nPrecision = %.20lf\n\n", delta);
  }
  init();

  build_antidictionary();

  if (OUTPUT_PROGRESS){
    fprintf(logfile, "Maximal queue size = %I64d\n", max_words_queue_size);
    fprintf(logfile, "The number of words = %d\n", words_num);
    fprintf(logfile, "The number of vertices = %d\n", total);
    fflush(logfile);
  }

  /*if (GUI_VERSION){
    if (use_long_mask) max_words_queue_size >>= 1;
    else max_words_queue_size >>= 2;
    printf("%I64d\n", max_words_queue_size);
    printf("%d\n", words_num);
    printf("%d\n", total);
  }*/

  build_dfa();

  if (OUTPUT_PROGRESS){
    fprintf(logfile, "\nDFA built!\n\n");
    fflush(logfile);
  }

  divide_scc_solve();

  fprintf(logfile, "Upper bound: %.13lf\n", grate);
  printf("%.13lf\n", grate);

  fflush(stdout);
  fclose(logfile);
  clear_data();
}

// Parses the file with parameters

void read_ini_file(){
  FILE * inif = fopen("params.ini", "r");
  char str[100], val[100];
  while (fscanf(inif, "%s = %s\n", &str, &val) > 0){
    if (!strcmp(str, "QUEUE_SIZE")){
      if (sscanf(val, "%d", &QUEUE_SIZE) > 0){
        auto_queue_size = false;
      }
    }
    if (!strcmp(str, "DFA_SIZE")){
      if (sscanf(val, "%d", &DFA_SIZE) > 0){
        auto_dfa_size = false;
      }
    }
    if (!strcmp(str, "MAX_ITER")){
      sscanf(val, "%d", &MAXITER);
    }
    if (!strcmp(str, "OUTPUT_PROGRESS")){
      OUTPUT_PROGRESS = (val[0]=='t' || val[0]=='T' || val[0]=='1');
    }
    if (!strcmp(str, "USE_ORACLE")){
      USE_ORACLE = (val[0]=='t' || val[0]=='T' || val[0]=='1');
    }
    if (!strcmp(str, "ADD_LOOPS")){
      ADD_LOOPS = (val[0]=='t' || val[0]=='T' || val[0]=='1');
    }
    if (!strcmp(str, "GEN_EDGE_LIST")){
      if (val[0]=='t' || val[0]=='T' || val[0]=='1') GEN_EDGE_LIST = 1;
      else if (val[0]=='f' || val[0]=='F' || val[0]=='0') GEN_EDGE_LIST = 0;
      else GEN_EDGE_LIST = -1;
    }
    if (!strcmp(str, "IMPROVE_TARJAN")){
      IMPROVE_TARJAN = (val[0]=='t' || val[0]=='T' || val[0]=='1');
    }
  }
}

int main(int argc, char * argv[]){
  read_ini_file();
  // File IO
  if (argc >= 2){
    freopen(argv[1], "r", stdin);
    if (argc >= 3){
      freopen(argv[2], "w", stdout);
    }
    while (scanf("%d%d", &ALPH_SIZE, &ROOT_LENGTH) > 0){
      if (ALPH_SIZE < 0) break;
      scanf("%d%d%d", &beta.num, &beta.denum, &beta.p);
      scanf("%d", &power_type);
      scanf("%lf", &delta);
      solve();
    }
  }
  // Or console mode
  else{
    printf("Hello! Working in console mode now...\n");
    while (1){
      printf("Alphabet size: ");
      scanf("%d", &ALPH_SIZE);
      printf("Max length of the root: ");
      scanf("%d", &ROOT_LENGTH);
      printf("Beta numerator: ");
      scanf("%d", &beta.num);
      printf("Beta denominator: ");
      scanf("%d", &beta.denum);
      printf("Is there a + ?: ");
      scanf("%d", &beta.p);
      printf("Type of abelian power (0 - weak, 1 - moderate, 2 - strong: ");
      scanf("%d", &power_type);
      printf("Delta (exponential form is acceptable): ");
      scanf("%lf", &delta);
      solve();
      printf("\nContinue? (0/1): ");
      int ch;
      scanf("%d", &ch);
      if (!ch) break;
    }
  }
  return 0;
}
